Ab initio molecular dynamics calculations of ion hydration free energies 



Kevin Leung* 
Surface and Interface Sciences Department, MS 1415 



o 
o 

1—5 

(N 
(N 

o 

S 

> 

a 

B 

o 
o 



> 
o 

o 

ON 

o 



X 



Susan B. Rempe 
Nanobiology Department, MS 0895 

O. Anatole von Lilienfeld 
Multtscale Dynamic Materials Modeling Department, MS 1322 
Sandia National Laboratories, Albuquerque, New Mexico 87185, USA 
(Dated: June 22, 2009) 

We apply ab initio molecular dynamics (AIMD) methods in conjunction with the thermodynamic 
integration (TI) or "A-path" technique to compute the intrinsic hydration free energies of Li^ , Cl~ , 
and Ag^ ions. Using the Perdew-Burke-Ernzerhof (PBE) functional, adapting methods developed 
for classical force field applications, and with consistent assumptions about surface potential (</)) 
contributions, we obtain absolute AIMD hydration free energies (AGhyd) within a few kcal/mol, 
or better than 4%, of Tissandier et al.'s [J. Phys Chem. A 102, 7787 (1998)] experimental values 
augmented with the SPC/E water model predictions. The sums of Li'^/Cl~ and Ag+/Cl~ AIMD 
AGhyd, which are not affected by surface potentials, are within 2.6% and 1.2 % of experimental 
values, respectively. We also report the free energy changes associated with the transition metal 
ion redox reaction Ag"*" + Ni^ — > Ag + Ni^^ in water. The predictions for this reaction suggest 
that existing estimates of AGhyd for unstable radiolysis intermediates such as Ni^ may need to be 
extensively revised. 
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I. INTRODUCTION 

Accurate predictions of hydration free energies of 
ions and molecules are crucial for modeling chemical 
and biochemical reactions in water and the adsorption 
of ionic species at water-material interfaces and inside 
nanoporesii State-of-the-art Density Functional Theory 
(DFT)-based ab initio molecular dynamics (AIMD) sim- 
ulations allow modeling the breaking and making of 
chemical bonds, as well as molecular polarizability. Di- 
rect use of AIMD to predict ion hydration free ener- 
gies, AGhyd, will have significant impact on computa- 
tional electrochemistry, biophysics, desalination, energy 
storage applications, corrosion studies, and geochemistry. 
AIMD simulations have already been extensively applied 
to study the hydration structure of ionSf^ii^i^i^ in many 
cases leading to more accurate predictions of the hydra- 
tion number than classical force field methods. At the 
same time, using hydration structure information and 
DFT and quantum chemistry calculations, the quasi- 
chemical method has been applied to predict highly accu- 
rate AGhyd for ions in water and biological binding sitesi^ 
In this manuscript, we generalize and apply AGhyd meth- 
ods developed for classical force fields to AIMD simula- 
tions. In some cases, our work can be related to "alchem- 
ical" potentials within the context of molecular grand- 
canonical ensemble DFT that allows variations of atomic 
numbers and electron numbers,^ 

Many of the techniques we use for predicting AIMD 
AGhyd have non-DFT precedents. In classical force field 



treatments of hydrated ions, AGhyd at infinite ion di- 
lution has been successfully computed^i^ i^"'^^ using the 
thermodynamic integration (TI) method j^^i^^ 



AGhyd — 



(1) 



or free energy perturbatiorii^ and closely related tech- 
niques. Here < A < 1 interpolates between the initial 
and final systems, H{X) is the Hamiltonian as A varies, 
the brackets denote equilibrium sampling with the Boltz- 
mann factor exp[— /3i7(A)], and (3 = I/UbT. For obvious 
reasons, the method is also called "A-path integration."— 
AGhyd is a state property, independent of the interpo- 
lation pathway. Force field parameters for ions are gen- 
erally fitted with a specific water model (e.g., SPC/E^^) 
to reproduce experimental AGhyd values. In simulations 
of monoatomic ions M with charge q, A is conveniently 
set to be proportional to q in Eq. [1] such that the ion is 
"charged up" linearly from M° to M'+ . 

Two critical theoretical advances have enabled di- 
rect comparisons of predicted AGhyd with tabulated 
data. (A) The long-range nature of coulomb interactions 
means a significant simulation cell size dependence arises 
when using Ewald summations. This dependence de- 
rives from the interactions of an ion with its images as 
well as with the neutralizing background in a charged 
simulation cell. To remove this dependence. Hummer, 
Pratt, and Garcia devised a monopole correction so ef- 
fective that even an 8-water simulation cell containing a 
Na'*" ion already yields AGhyd well converged with sys- 
tem size.^'^ (B) Comparison with experiments effectively 
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entails bringing an ion from vacuum at infinity into the 
bulk liquid water region. A surface potential, (j), material- 
izes at the liquid-vapor interface, leading to a shift in the 
ion free energy q4) in the aqueous phasCfl^J^d- Account- 
ing for the surface potential, the calculated absolute ion 
hydration free energy, which may not be measurable jiZ, 
becomes 

AGtot = AGswaid + q{<j)d + (2) 

Here AGswaid is the hydration free energy computed us- 
ing standard Ewald summation which assumes a zero av- 
erage electrostatic potential inside the simulation cell.^" 
(j)d and (j)q are the dipolar and quadrupolar (or "spherical 
second moment") contributions to the surface potential 
4>. Some reported experimental data have subtracted the 
effect of this potentiaP^ while others have uot^ 

The rapid convergence of AGhyd with simulation cell 
size (A) significantly facilitates the application of this 
AGhyd formalism to computationally costly DFT-based 
AIMD simulations. Special attention should be paid 
to the surface potential contribution (B) in AIMD set- 
tings. Unlike classical models for water, (f) — 4>d + <i>q 
has not yet been predicted for AIMD water (e.g., com- 
puted with a generalized-gradient approximated (GGA) 
Kohn-Sham functional such as Perdew-Burke-Ernzerhof 
(PBE)^'^). Such a calculation would entail a large sim- 
ulation cell depicting the interface and long sampling 
trajectories. Furthermore, as the liquid water density 
affects 0^ jiiii^iiSi^l the effectiveness of such a calcula- 
tion may further be limited by the fact that bulk GGA 
water may not exhibit 1.0 g/cm^ density.— Although 
4id and (t)q are not independent — they require a com- 
mon choice of molecular center, typically taken to be 
the oxygen atom of water molecules — the quantity (j)q 
has recently been computed for PBE water using maxi- 
mally localized Wannier functionSf^^ This piece of infor- 
mation is important for DFT-based calculations because 
AGEwaid itself is an ambiguous quantity whose value de- 
pends on whether the pseudopotential contains core elec- 
trons, while AGswaid + 90g is independent of such DFT 
details. We therefore redefine 

AGhyd = AGEwaid + q4>q- (3) 

To further enable comparison with experimental data 
in Ref. [2l|, which contain no surface potential contribu- 
tions, we add q(f>q~-19.7q kcal/mol, the quadrupole mo- 
ment value for SPC/E water at 1.00 g/cc density when 
the oxygen site is chosen as the molecular center. This 
is appropriate because AGhyd for various ions have been 
fitted to Ref. 21 using the SPC/E water model^o or the 
very similar SPG model.— ^S, In effect, we are comparing 
AIMD AGhyd with SPC/E calculations fitted to the data 
of Ref. dJ. For the data tabulated in Ref. 2^, which 
contain the surface potential term q{(l)d + 4>q)i sub- 
tract (70d = 4.8q kcal/mol estimated using SPG/E water 
model-based water-vapor interface molecular dynamics 
calculations Although an investigation of 0^ predicted 



with different methods is not the focus of this work, accu- 
rate DFT methods and accurate force fields should yield 
similar, reliable (fjd. Even if there exists a 50% uncer- 
tainty in this SPG/E (pd estimate, AGhyd + <l4'd in wa- 
ter will be affected by only ~ 2.4|g| kcal/mol. Indeed, 
the much used SPG and the TIP4P water models yield 
(^ci=5.5 and 7.1 kcal/mol/|e|, respectively t^^'^^i^^ which 
are slightly different from the SPG/E (t>d. The discrep- 
ancies among these models can be taken as a measure 
of the systematic uncertainty associated with our (j)d as- 
signments. 

Finally, experimental data for moving ions from vac- 
uum into aqueous solution are referenced to their respec- 
tive standard states, i.e., gas phase ions at 1.0 atm. pres- 
sure and hydrated ions at 1.0 M concentration. To be 
consistent with the infinite dilution limit AGhyd pre- 
dicted in this work, G*^°^ — 1.9 kcal/mol is further sub- 
tracted from tabulated AGxiss for all ions regardless 
of their charges to account for the volume change in- 
cluded in the experimental data. Due to a sign problem,— 
2G^°^ kcal/mol needs to be subtracted from AGMarcus for 
this purpose. 

To summarize, we compare our AIMD AGhyd (Eq- El 
with AGMarcus + 90g^'"^^ — 2G("^ kcal/mol and AGxiss — 
qcj)^^'"^^ —C^^^ kcal/mol, where AGMarcus and AGxiss are 
the values listed in Refs. |2ll andlU, respectively. 

Note that the proton is often used as a reference for hy- 
dration free energies.— Referencing the predicted AGhyd 
of ions with that of H"*" computed in the same way cir- 
cumvents the need to estimate 0. In AIMD settings, 
however, an excess proton can migrate from one II2O to 
another. Therefore we have not yet attempted to com- 
pute this proton AGhyd- 

For test cases, we consider Li+ and Gl^. The Li+ 
ion hydration structure and hydration free energies have 
been extensively studied using AIMD and quasi-chemical 
methods, respectivelyjS, Gomputing the AGhyd of Gl~ 
further allows us to predict the summed AGhyd of the 
monovalent Li+/Gl~ pair, where the surface potential 
terms cancel and the result contains less systematic un- 
certainty. We show that this summed value is at worst 
within 2.6% of experimental resultsi^ii^ 

We also study the change in hydration free energies 
associated with 

Ag+ + Ni+ ^ Ag + , (4) 
and the corresponding electrochemical half-cell reactions, 
Ag Ag+; and (5) 

Ni+ m^+. (6) 

These reactions are pertinent not only to elementary elec- 
trochemical processes, but also to the initial stages of 
nano alloy synthesis by radiolysis^ii^ 7 irradiation of 
mixed electrolytic aqueous solutions releases secondary 
electrons that reduce the metal ions to atoms or lower ox- 
idation state ions. These reduced species readily coalesce 



to form clusters. In the case of a mixed Ag(I) /Ni(II) so- 
lution, the exothermicity of Eq.|4]will determine whether 
reduced Ni species are readily re-oxidized by Ag"*" in the 
solution — a side reaction that hinders nano-alloy cluster 
formation. AIMD is an attractive route to estimate the 
redox free energies associated with Ni(I) species, which 
exhibit short lifetimes and are difficult to probe experi- 
mentally. 

Apart from the ability to compare AIMD AGhyd with 
quasi-chemical theoryS*^ and potentially extend DFT- 
based absolute hydration free energy calculations to in- 
homogeneous media, this work is important due to its 
close relationship to recent theoretical advances. One 
is the alchemical A-path integration technique recently 
formulated within a DFT/AIMD-based molecular grand 
canonical ensemble scheme^^ which accounts for changes 
in pseudopotentials as well as the number of electrons. 
As long as the pseudopotential replaces all core electrons 
in the ion, AGhyd TI calculations are very similar within 
AIMD and the SPC/E model treatments of water. More 
complex treatments are required, however, when ion in- 
sertion into the solvent involves not only changes in the 
ionic pseudopotential, but also injection of electrons.^ 
This alchemical path technique has been applied to quan- 
tum mechanics/molecular mechanics (QM/MM) simu- 
lations of electron transfer reactions of aqueous metal 
complexes (Fe(II/III) and Ru(II/III)).^ Our work is 
even more closely related to purely AIMD-based com- 
putational electrochemistryi^ Here the electron transfer 
processes are similar to those in Ref. [sj, but all wa- 
ter molecules are treated with DFT methods, and the 
long-range electrostatics are fundamentally different from 
those in QM/MM calculations. Our computational ap- 
proach treats the ionization potential and the ion hy- 
dration free energy contributions to the redox potential 
separately. While it is based on and derives its rigor from 
theories well established with classical force field hydra- 
tion treatments (e.g., Eq.[3]), our thermodynamic method 
has not been extended to estimate the fluctuating gaps 
that are necessary for calculating reaction rates via the 
Marcus theoryi^^ 



II. METHOD 

A. VASP calculations 

We apply the Vienna atomistic simulation package 
(VASP)^ version 4.6 with a modified pot.F,i^ the PBE 
exchange correlation functional, projected-augmented 
wave (PAW) pseudopotentials'^^'^^ (PP) with only va- 
lence electrons for Li, CI, H, and O atoms, and Ag and Ni 
PPs that include pseudovalent 4p and 3p electrons. Two 
protocols to generate VASP AIMD trajectories for Li+ 
solvated in water are applied. For the ion plus 32-water 
simulations, we use a cell size of 9.855 A corresponding to 
a water density of 1.0 g/cc, a 0.25 fs time step, an energy 
cutoff of 400 eV, and a Born-Oppenheimer convergence of 
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FIG. 1: The binding energies and optimized distances be- 
tween a H2O molecule and VASP PBE pseudopotentials glob- 
ally scaled by a factor of < g < 1. (a) & (c): Li+; (b)&(d): 
Na^. The pseudopotentials have no core electrons. Dashed 
lines are cubic spline fits. Na^ is meant as a counter example 
to Li"*" for gas phase behavior; its behavior in water will not 
be the focus of this work. 



10 ^ eV at each time step. For 64- water simulations, the 
corresponding parameters are 12.417 A (1.0 g/cc), 0.5 fs, 
500 eV, and 10~^ eV, respectively. These settings limit 
the temperature drifts to 1 and 0.5 K/ps, respectively. 
The trajectory length for each value of q is at least 40 ps 
in 2-point TI calculations and at least 30 ps for 6-point 
TI. Initial configurations are pre-equilibrated using the 
SPC/E water model and ion force fields^" with charges 
scaled to the net charge of the corresponding AIMD 
simulation cells. A Nose thermostat is applied, setting 
T=400 K, which is needed for the PBE functional to de- 
scribe experimental liquid water at room temperaturei^ 
The deuterium mass is adopted for all protons to allow a 
larger time step, although the H mass is assumed when- 
ever water density is reported. Ag"*" and Ni^+ simulations 
are performed at 0.99 g/cc water density using 32-II2O 
simulation cells. 



B. Visualizing Electronic Isosurfaces 

Electronic isosurfaces and integrated changes in elec- 
tron density, A(a;) = J dydz[p{x,y, z)n ~ pix,y,z)c] as 
functions of spatial coordinate x, are also computed and 
depicted for Li''"'" in water for various values of q. The de- 
picted geometries are snapshots taken at the end of the 
32-water PBE simulations. These results are obtained 
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FIG. 2: {dH{q)/dq)q for Li^"*" as q varies. The bare ion 
contributions, Ewald corrections, and electrostatic potential 
shift due to the quadrupole moment have been subtracted. 
Crosses: 32 H2O, 1.00 g/cc; circles: 64 H[20, 1.00 g/cc; trian- 
gles, same as crosses but are for SPC/E water. The dashed 
lines are cubic least-squared fits to the crosses and triangles. 
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FIG. 3: (a): Integrated changes in electron density, A(x) = 
J dydz[p{x,y, z)n — p{x,y, z)c\, as a function of spatial coordi- 
nate X for the various values of q. p„ and pc are the densities 
for the neutral and the charged systems, respectively. All 
charged species, Li'"*" have been shifted to a; = 0. Symbols 
correspond to actual grid-points, the continuous lines are cu- 
bic interpolations, (b)-(d): Isosurface plots of the electron 
density difference, p(x, y, z)n — p{x, y, z)c (iso-value = ± 0.01 
a.u., white < 0, blue > 0), for q =0.1, 0.6, and 1.0. Periodic 
boundary conditions apply; the prominent, 8 blue spheres rep- 
resent the (periodically replicated) changes in Li'^ densities, 
and some changes in water dipole moments are apparent too. 
See Sec. Ill Bi for technical details. (Panels (b)-(d) are located 
at the end of the document.) 



using the code CPMB^P the PBE functional,— pseudopo- 
tentials from Ref. [M and a cutoff of 100 Ry (1361 eV). 
Pc refers to the electron density obtained by minimizing 
the energy within the indicated charge. As with VASP, 
CPMD uses an opposite background charge to neutralize 
the system within the periodically replicated simulation 
cells, pn corresponds to the density of the same geom- 
etry but with the charged species replaced by a neutral 
He atom. 



C. Li+ thermodynamic integration 

To implement Eq. [T] for Li"*", we generate integrand 
values at different q values according to two different 
integration formulas: a two-point Gaussian quadrature 
and a six-point trapezoidal rule. To that end, AIMD 
trajectories apply a Li"*" pseudopotential (which contains 
no core electrons) globally scaled by Gaussian quadra- 
ture values g=0.211325 and 0.788675. This procedure 
is analogous to the scaling of the ionic charges in clas- 
sical force field molecular dynamics calculations of hy- 
dration free energieSf^S In addition, 9=0.1, 0.4, 0.6, and 
1.0 are considered. Using these 6 points, a cubic least- 
squared fit is applied to extrapolate the integrand value 
to (7=04^ These steps yield 6 almost evenly spaced in- 
tegration points needed to implement a trapezoidal rule 
integration. 



Figure [T^ shows that the scaled VASP Li+ pseudopo- 
tential behaves to some extent like a classical force field 
Li*+; its binding energy with one H2O molecule scales 
roughly linearly with q except at very small q. The op- 
timal Li-Owator distance also shrinks smoothly with de- 
creasing q (Fig. [T};). In contrast. Fig. [T|3 shows that the 
scaled VASP PBE Na+ exhibits water binding energies 
that deviate more strongly from linearity. Furthermore, 
the optimal q-scaled Na+-0H2 distance sharply decreases 
to 0.87 A at g w 0.29, which suggests the formation of an 
anomalous covalent bond beyond q < 0.29 (Fig.[TJi). For 
efficient AIMD AGhyd simulations, a pathway should be 
chosen such that at the selected simulation points, elec- 
tron transfer or unphysical chemical bonding between the 
scaled pseudopotential and H2O is avoided. 

The AIMD trajectory is sampled every 0.1 ps. At such 
intervals, we use a finite difference method to compute 
dH{q)/dq = [H{q + Aq/2) - H{q - Aq/2)]/Aq at fixed 
atomic configurations. Here H{q) is the total potential 
energy of the simulation cell predicted using VASP. When 
taking finite derivatives, Aq values of 0.025 and 0.050 
yield Li+ hydration free energies that agree to within 
0.5 kcal/mol. Evaluating {dH{q)/dq)q using 400 eV and 
500 eV cutoffs lead to indistinguishable results. 

The derivative is corrected for finite size effects by 
adding the Ewald correction to the energy aq^/2L at 
each q, where a is the Madelung constant, to the 
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Li+-plus-water VASP energies (issue "(A)" discussed in 
the introduction^). The quadrupole moment correc- 
tion q(l)q is hnearly dependent on q and has been esti- 
mated in Ref. J9,. With the sUghtly smaller simulation 
cell used in this work, the (f)q corrections are predicted 
to be 3.85 and 3.81 eV for 1.00 and 0.99 g/cc water 
density4^ Unlike classical force field calculations, the iso- 
lated ion Li''+ carries a non-zero energy. Thus we sub- 
tract (d-ffbarc ion 

{q)/dq)q from Eq. [H 
Unless otherwise noted, the Li+ thermodynamic inte- 
gration protocol (e.g., the sampling interval, subtraction 
of bare ion energies) is applied to all other ions. 

D. Cl^ thermodynamic integration 

AGhyd for Cl^ requires a different TI procedure. Un- 
like the Li+ PP without explicit Is electrons, scaling 
the VASP Cl~ PP to zero also involves removing 8 elec- 
trons. While it is possible to alchemically perturb Cl~ 
to Ar, this TI route is not directly applicable for multi- 
atom anions. Instead, we first use TI to "grow" a non- 
polarizable classical force field (FFl^" Cl~ with a nega- 
tive point charge and a Lennard- Jones interactiorti^ with 
the oxygen sites of PBE water. This can be regarded as 
a QM/MM simulation, but with the solvent (not solute) 
treated quantum mechanically. Then we use a one-step 
free energy perturbation (FEP) procedure, 

/3[AG(PBE) - AG(FF)] 
= -log(exp[-/3(i/(PBE)-ff(FF)])FF, (7) 

to estimate the PBE Cl^ AGhyd- As long as the hy- 
dration structures of the classical and PBE ion in PBE 
water are similar, this method can be generally and accu- 
rately applied to multi-atom anions or cations, as well as 
PP's like the VASP PAW PBE Na+ whose interaction with 
water exhibits anomalies when the PP is scaled continu- 
ously to zero (Fig. [T]) . If there are partial positive point 
charges in the classical force field, however, the DFT va- 
lence electrons may collapse onto those atomic sites, and 
pseudopotentials that repel electrons may be needed to 
prevent such a collapse. 

E. Ag+ and Ni^^ thermodynamic integration 

The VASP PBE pseudopotentials used for Ag and Ni 
contain 11 and 16 electrons, respectively. When the 
number of electrons in 32-water simulation cell is fixed 
at (32 X 8 11 - q) and (32 x 8 16 - q) in AIMD 
trajectories, our maximally localized Wannier function 
analyses^ reveal that (11 — q) and (16 — q) electrons re- 
main localized on Ag and Ni, respectively. This indicates 
that Ag'^+ and Ni'^+ species exhibit no tendency to eject 
excess electrons into water )^ and the partially charged 
ions are preserved within a A-path that vary the total 
number of electrons in the system. Hence we simply use 
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TABLE I: Li+ hydration free energies using different compu- 
tational protocols. H2O densities and AGhyd are in units of 
g/cc and kcal/mol, respectively. "Ref. [U; ''Ref. [S^. Exper- 
imental values adjusted for surface potentials and standard 
state contributions are marked with a dagger (see text). 



the number of electrons as the order parameter. A, analo- 
gous to Refs.HjIH, andim is simply computed by 
adding and subtracting 0.025 electrons to the simulation 
cell and performing a finite difference. The exceptions 
are Ag"*" (where we compute the difference between Ag"*" 
and Ag0-95+); Ni+ (Ni+ and Nii-05+); and Ni2+ 
and Ni2+). As we subtract the bare ion contribution at 
each q, the expression ((^^^^) — i^^^^^i^-i^^iil ) should re- 
flect purely solvent-induced effects. 



For Ag, spin-polarized PBE calculations are adequate. 
In contrast, spin-polarized PBE-based AIMD simulations 
of Ni'"*" in water underestimate the gap between the high- 
est occupied (HOMO) and lowest unoccupied (LUMO) 
molecular orbitals. This occurs because PBE severely 
underestimates exchange interactions in the localized 3c? 
orbitals, leading to near degeneracies in intermediate- 
q Ni''"'" d-shell orbitals and slow numerical convergence 
of the electronic structure at each Born-Oppenheimer 
AIMD time step. We have therefore applied the DFT-I-U 
technique''^ to the Ni 3d orbitals to generate AIMD tra- 
jectories with which we evaluate Eq.[T]using only the PBE 
functional. Originally devised for solid state applications, 
DFT-I-U has recently been adapted for molecular systems 
and even used in AIMD settings.'*^:'*® U is set at 4.0 eV 
to yield a 15.7 eV gas phase Ni^"*" binding energy in a 
Ni^+(H20)6 cluster. This is the value predicted using the 
B3LYP hybrid functional*^ and a 6-311-|-G(d,p) basis.^S 
Using DFT+U generated geometries for PBE AGhyd is 
justified because, in the gas phase, the PBE functional 
and DFT-I-U predict optimized Ni^"'"(H20)6 geometries 
which are nearly identical. 



III. RESULTS 



A. Li^ hydration free energy 

Figure O plots {dH (q) / dq) q as q varies after subtract- 
ing contributions from Ewald images;^ the quadrupole 
or spherical second moment contribution qcjjq^r- and the 
energies of the bare Li^+. {dH{q)/dq)q computed using 
32- and 64-H2O simulation cells at 1.00 g/cc H2O den- 
sity are in good agreement at q = 0.21 and q — 0.79. 
Using a 2-point Gaussian quadrature, AGhyd for the two 
cells integrate to -128.6 and -126.7 kcal/mol, respectively 
(TablelU. Splitting the data into four segments, the stan- 
dard deviations in these AGhyd are found to be 1.1 and 
0.5 kcal/mol, respectivelyj^ Thus the two cell sizes ex- 
hibit AGhyd approximately within numerical uncertain- 
ties of each other, showing that the finite system size 
effect is small for AIMD after applying the Ewald correc- 
tion, as is the case with classical force field simulations]^ 
A dielectric continuum estimate would suggest that, after 
adding the leading order (1/i) Ewald correction, the 32- 
water simulation cell is already converged to the infinite 
dilution limit to within 1 kcal/moli^ 

For illustrative purposes, we also display in Figure [3] 
the Li-ion growth-induced changes of the total electron 
density integrated over the x- and j/-coordinates. From 
inspection of this change arising from the presence of 
the increasingly charged ion one can conclude that, as 
expected, the attraction of electrons toward the ion in- 
creases as the charge approaches +1.0. The isosurfacc 
plots support a similar conclusion. For small values of 
q, changes in density occur throughout the system. As q 
approaches its final value, however, the drastic increase 
in electronic density at the ion position due to increas- 
ingly polarized water (Fig. [3^) is hidden behind the large 
sphere of depleted density. This large sphere comes about 
because we have subtracted the electron density of a neu- 
tral helium atom from that of the Li"*" pseudopotential. 

Figured] depicts the pair correlation functions g{r) be- 
tween Li''"'" and the O and H sites in H2O. Recall that 
the entire VASP PBE Li"*" pseudopotential, including the 
long-range coulomb and the short-range Pauli-exclusion 
contributions, is scaled with q. Hence, at small g, the 
most probable Li'"'"-Owator distance is much reduced from 
the 5 = 1 case. Nevertheless, we have verified that neg- 
ligible electron density resides near the Li'"*" nuclei, indi- 
cating that Li''"'" does behave like a partially charged ion 
in water. The insets depict the instantaneous hydration 
numbers A^^,, computed at each time step by integrat- 
ing each 5Li-o('') to its first minimum. For q = 0.21, 
Nyj averages only to 1.5 and experiences rapid temporal 
fluctuations. Despite this, 5Li-o(?') still exhibits a high 
peak value because the scaled Li''"^ has such a small ra- 
dius. At g = 0.79, — 3.5, approaching the = 4 
AIMD value reported for Li"'".^ 

Figure [S] depicts the logarithm of the distributions of 
instantaneous hydration numbers for Li'^'^"'" and Li"*". In 
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FIG. 4: Pair correlation functions g{r) between Li'"*" and the 
O (solid line) and H (dashed line) sites of H2O molecules, (a) 
g=0.21; (b) 5=0.79. The instantaneous hydration numbers 
are depicted in the insets. 



conjunction with low order m {d™ H (q) / dq™) derivatives, 
hydration number distributions at the TI end-points can 
in princple be used to predict the hydration free energy 
using a single AIMD trajectory a,t q = 01 q = 1»^ Since 
we have avoided q = and the finite differences applied 
in our implementation may not be accurate for m > I, we 
have not attempted to estimate AGhyd with high order 
derivatives, but have used 2 or 6 q values to evaluate 
AGhyd- Note that, using the quasi-chemical theoretical 
framework, hydration number distributions of a solute 
can be used directly to estimate hydration free energies;^, 
as demonstrated in recent works-^'^' -'* Furthermore, such 
distributions are of intrinsic interest and can lend useful 
comparison with those predicted using classical force field 
simulations. See also Ref. [s^ for other methods devised 
to reduce the number of q-value integrands needed to 
perform TI calculations. 

We next investigate the accuracy of the 2-point TI 
quadrature by further sampling {dH(q)/dq)q at 17=0.1, 
0.4, 0.6, 1.0 in addition to 0.21 and 0.79 in a simula- 
tion cell. This denser grid allows an approximate 6- 
point trapezoidal rule integration after we extrapolate 
{dH{q)/dq)q to q=0.0. Figure [2] shows that {dH{q)/dq)q 
is almost linear for a large, intermediate q range except 
near g=0 and <7=1. This is in qualitative agreement 
with SPC/E model prediction o^'^^ which we also com- 
pute for a 32-water simulation cell and depict in Fig. [2l 
The deviation from linearity at g = is well-reproduced 
with a cubic fit for both AIMD and SPC/E {dH{q)/dq)q. 
Table U confirms that the 2-point and 6-point formulas 
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FIG. 5; Logarithm of the probabihty (Pn) of instantaneous 
hydration numbers (n) multiphed by thermal energy, in units 
of kcal/mol. (a) Li''"'"; (b) Ag'"*". Squares and dashed lines: 
q = 0.2; circles and solid lines; g = 1.0. n is determined 
by counting all water oxygen atoms within 2.08, 2.75, 2.90, 
and 2.92 A of the four ions, respectively. These distances are 
determined by locating ithe first minimum in the ion-water 

air). 
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FIG. 6: (a),(b) Pair correlation functions g{r) between classi- 
cal force field Cl'~ and the O (solid line) and H (dashed line) 
sites of PBE H2O molecules, (a) g=0.21; (b) (7=0.79. (c) 
{dH{q)/dq)q for classical force field CI'" as q varies. Crosses 
and triangles are for AIMD and classical force field treatments 
of water in 32-H2O simulation cells. The bare ion contribu- 
tions, Ewald corrections, and electrostatic potential shift due 
to the quadrupole moment have been subtracted. The dashed 
lines are cubic least-squared fits. 



yield AGhyd within 0.3 kcal/mol of each other — well 
within the numerical uncertainties of the simulations. 
Henceforth we will report the 6-point value of AGhyd 
= -128.3 ± 0.9 kcal/mol for Li+. 

This success of the 2-point formula appears however 
somewhat fortuitous. One would not a priori expect this 
quadrature to be accurate for Li"^ because of the large 
changes in effective Li''+ radius (Fig. 3] below). The 
classic Born hydration free energy formula, based on a 
dielectric continuum description of the solvent, predicts 
AGboiii oc g^/(2a)(l — 1/e) at a fixed ionic radius a. It is 
quadratically dependent on q when a is held constant. In 
non-polarizable classical force field AGhyd simulations, 
the Lennard-Jones radius of the ion is also held fixed 
while the ionic charge varies. The constant radius thus 
seems crucial to the accuracy of the 2-point Gaussian 
quadrature, which is exact only if {dH{q)/dq)q is linear 
in q. Despite this, the 2-point formula will be shown to be 
accurate for the AIMD AGhyd associated with Li+ , Ag+ , 
and Ni"*" — > Ni^"*" considered in this work. It appears less 
accurate for Cl^, unlike SPC/E-based Cl^ AGhyd calcu- 
lations. The fact that the radius of Li''"'" (and to some 
extent, other ions) changes with q in our DFT calcula- 
tions also explains the discrepancy between AIMD and 
SPC/E (dH{q)/dq)g^a values. 



To compare AIMD predictions with experimental data, 
AGMarcus + <Z</'g^'^^^ ~ 2C^'^^ kcal/mol is found to be - 
137.0 kcal/mol,2i while AGTiss-^^^^^-C^"^ kcal/mol 
=-133.2 kcal/mol22 (TablelJ). These values are similar to 
the SPC/E AGhyd for Li"^, and are 8.7 and 4.9 kcal/mol 
higher than the 6-point AIMD prediction for a 32 H2O 
simulation cell, respectively. The discrepancies with 
AIMD predictions may be due to numerical noise, PBE 
functional inaccuracies, or systematic uncertainties aris- 
ing from the treatment of \e\(j)- Indeed, the discrepancy 
between SPC/E-augmented experimental values listed by 
Marcus2i and Tissandier et al^ can also be taken as a 
measure of surface potential-related systematic ambigu- 
ity. This issue will be interrogated in the next subsection 
when we consider the anion Cl"~. 

An optimal study of hydration free energy would in- 
clude also the changes in water density due to the pres- 
ence of salt cations and anions or water confinement in- 
side nanopores. We have therefore examined the effects 
of reducing the water density to 0.97 g/cc. This small 
reduction in water density corresponds to the activity 
of water at 0.1 M ion concentration, which is the typi- 
cal concentration of K"*" ions in the cytoplasm of skeletal 
muscle cells and the typical concentration of Na"*" and 
CI"' ions outside cellsj^. Table U shows that the small 
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TABLE 11: CI" hydration free energies. The asterisk denotes 
AIMD AGhyd adjusted for finite simulation cell size and pack- 
ing effects (see text). Also listed are AGhyd for Li+ plus CI". 
The SPC/E results for CI" and Li+/Cl" contain the packing 
correction. H2O densities and AGhyd are in units of g/cc and 
kcal/mol, respectively. "Ref. [21I: ''Ref. [2^. Experimental val- 
ues adjusted for surface potentials are depicted with a dagger; 
see text for details. 

effect on AGhyd due to water density changes is within 
the numerical uncertainty. This weak dependence is con- 
sistent with quasi-chemical theory analysis^ where con- 
tributions to AGhyd are separated into inner hydration 
shell and outer shell contributions. In the "cluster" im- 
plementation of the theory;^ the former can be deter- 
mined from gas phase cluster calculations scaled by water 
density, while the latter depends on the water dielectric 
constant, which is relatively independent of H2O density. 
As pointed out by Varma and Rempe^^ since the depen- 
dence of free energies on water concentration is logarith- 
mic, large changes in water density are required before 
there is an effect on AGhyd- 



B. CI hydration free energy 

Figures [5^ and b depict the g{r) between the classical 
force field CI'" (henceforth FF-CF") and the oxygen and 
proton sites of H2O molecules at two q values. At g=0.21 
(or even q=OA), FF-Cl'" is predominantly a hydropho- 
bic sphere that excludes both O and H from its vicinity. 
Due to the sheer size of the Lennard-Jones sphere that 
represents CF", this solute is seen to substantially dis- 
rupt the water structure around it in the 32-H2O simula- 
tion cell. Thus, in panel (b), the Cl-0 g{r) has dropped 
below 0.5 density units at r ~ 5 A — unlike the case 
for Li''+ at small q (Fig. [3^). At g=0.79, the ion forms 
hydrogen bonds with water; its gci-ii{r) exhibits a peak 
at r = 2.2A. At q=l (not shown), we obtain a FF-Cl''" 
hydration number of N^^bA, in good agreement with 
fuU AIMD simulations of PBE Cr in FEE water 



Figure[5J: depicts the variation of {dH{q) /dq)q FF-CP^ 
in PBE water as q varies. To obtain AGhyd for the PBE 
CI" ion, we further apply Eq. [7| to configurations sam- 
pled 0.1 ps apart along the AIMD trajectory. The dif- 
ferences between the instantaneous potential energies for 
FF-Cl" and PBE CI" are found to be almost constant 
with an estimated standard deviation of 0.15 kcal/mol. 
This indicates that FF-Cl" is an excellent reference for 
the PBE CI" . After a cubic polynomial extrapolation to 
q=0 and applying a 6-point integration formula, AGhyd 
for the PBE CI" integrates to -76.6±0.4 kcal/mol (Ta- 
ble |T]). A 2-point Gaussian quadrature formula yields 
-79.0±0.8 kcal/mol. As the latter is only exact for linear 
{dH{q)/dq)q, deviation from linearity in Fig. [S] indicates 
that a denser grid may be needed despite the constant 
radius of the FF-Cl sphere. This slight non-linearity is 
apparently due to water polarizability; corresponding 6- 
point and 2-point SPC/E calculations in 32- water simula- 
tion cells yield indistinguishable results. As {dH{q)/dq)q 
is well-fitted to a cubic polynomial in q and the trape- 
zoidal integration rule is accurate for cubic polynomials, 
however, Fig.[5J; strongly suggests that an integration for- 
mula higher order than the trapezoidal rule is not needed. 
Henceforth we report the 6-point TI value. 

Two post-processing corrections for AGhyd, unneces- 
sary for Li"*", need to be included here. (1) While the Li"*" 
PP is globally shrunk to zero, at ^=0 FF-Cl'^ remains 
a Lennard-Jones sphere that displaces water. This gives 
rise to an entropic or "packing" penalty; the contribution 
is estimated to be 4.0 kcal/mol using SPC/E water model 
simulations. (2) Simulation cell size effects are more sig- 
nificant for Cl^ than for Li+ , presumably because of the 
size of the Cl"^" sphere at small q (Fig.[S^). When we per- 
form purely classical force field simulations of a CI" ion 
in SPC/E water, we find that a 32-II2O simulation cell 
overestimates AGhyd by 3.3 kcal/mol compared to a 255- 
H2O cell. This discrepancy is much larger than the nu- 
merical uncertainty. In contrast, these two cell sizes yield 
Li"*" AGhyd that are within 1 kcal/mol. The simulation 
cell size dependence has been estimated using a dielectric 
continuum approach in Ref. [^. Assuming AIMD exhibits 
CI" packing penalty and simulation cell size dependence 
similar to classical force field MD, we add a 7.3 kcal/mol 
correction to the AIMD result. The corrected AIMD CI" 
AGhyd is listed in Table HIl It is within 0.4 kcal/mol of 

AGxiss — 'Z'/'d^^^^ ~ kcal/mol, and overestimates the 

magnitude of AGMarcus + q<P^'"^^ — 2G(°^ kcal/mol by 
4.0 kcal/mol. 

Adding AGhyd of oppositely charged monovalent ions 
eliminates the systematic uncertainty due to surface po- 
tential contributions. The combined AGhyd for Li+ and 
CI" are within 4.7 and 5.3 kcal/mol of experimental data 
quoted in Table |TT] respectivelyj ^^'^^ they underestimate 
those values only by about 2.3 and 2.6%. This sum, de- 
rived from MarcusSi and Tissandier et al.^ are within 
0.6 kcal/mol of each other, unlike in the cases of the iso- 
lated Li^ and CI" ions where the two adjusted experi- 




FIG. 7: Pair correlation functions g{r) between Ag'^^ and the 
O (solid line) and H (dashed line) sites of H2O molecules, (a) 
9=0.21; (b) g=1.00. The instantaneous hydration numbers 
are depicted in the insets. 



FIG. 8: {dH{q)/dq)g for Ag«+ and Ni«+ as q varies. The bare 
ion contributions, Ewald corrections, and global shift in the 
electrostatic potential due to the quadrupole moment have 
been accounted for. The dashed lines are cubic least-squared 
fits. 



mental data sets disagree by 3.8 and 4.4 kcal/mol, respec- 
tively. This suggests that the rather large, 8.7 kcal/mol 
discrepancy between AIMD AGhyd and Marcus' data for 
Li+ is partly due to the assignment of the SPC/E (j)q 
contribution to the surface potential. In contrast, Tis- 
sandier et aZ.'s data for the isolated ions are in substan- 
tially better agreement with AIMD AGhyd for both ions, 
suggesting that augmenting AGxiss with SPC/E is a 
reasonable approximation. 



C. Ag ^ Ag+ 

In Fig. [71 Ag-Owator and Ag-Hwatcr g{r) are depicted 
for two selected values of q. Unlike Li, the Ag atomic core 
is not scaled with q, and Pauli repulsion ensures that no 
water molecule penetrates the Ag core region. Thus the 
g{r) is not sharply structured at small q, and Ag^+ re- 
sembles a hydrophobic sphere as q decreases. For both q 
points, II2O in the first hydration shells are highly labile; 
see the insets. The Ag+-H20 g{r) (Fig. [Tja) yields a first 
shell hydration number N^—SA. The instantaneous hy- 
dration number distribution is depicted in Fig. [5l This 
is qualitatively similar to the Nw—4:-0 computed us- 
ing AIMD and another exchange correlation functional. 
Both these AIMD iV^, values are in good agreement with 
experimentS] ^°'^^ In contrast, a recent classical force field 
model with parameters fitted to quantum chemistry cal- 
culations has reported N^j = 6.®^ With the corrections 
(A)-(B) discussed earlier, a 6-point trapezoidal rule inte- 



gration, and a 1.6 kcal/mol packing correction estimated 
using classical force field simulations, we obtain AGhyd 
119.8±0.4 kcal/mol. This magnitude is 6.4 kcal/mol 

smaller than AGMarcus + q4'T''^^^ ~ 2G'-°'' kcal/mol (Ta- 
ble Unj-^i The sum of AIMD Ag+ and CI" AGhyd, how- 
ever, underestimates the experimental datap^ by only 
2.4 kcal/mol, or by 1.2 %. 

D. Ag+ + Ni+ ^ Ag + Ni2+ 

The details of Ni'+ hydration will be described 
elsewherci^ Here we focus on the change in AGhyd as 
Ni+ loses an electron. We use the PBE functional to com- 
pute {dH{q)/dq)q at 0.1 ps intervals along the DFT-fU 
AIMD trajectory with U=A eV. Figure [Hh shows that 
{dH{q)/dq)q is fairly linear as q varies. With a 6-point 
trapezoidal rule integration, Eq. [T] yields a change in 
AGhyd of -365.5±1.0 kcal/mol. A 2-point integration 
predicts a similar -363.4±2.4 kcal/mol. Unlike the cal- 
culations for Li"*" and Ag"*" , this system benefits from the 
fact that at "A"=(g - 1)=0, Ni+ is still highly charged, 
and larger statistical uncertainty at small q is avoided. 
Nevertheless, due to the slower water dynamics around 
the more highly charged Ni'^"'" ion, sampling correlation 
times may be longer and our error bars for Ni^+ may be 
underestimated. 

The electrochemical half cell reaction free energy con- 
sists of the change in AGhyd plus the ionization poten- 



10 



ion 






quadrature 


^ ily <_l 




32 


0.99 


2-pt 


-121.3 




32 


0.99 


6-pt 


-121.4 




32* 


0.99 


6-pt 


-119.8 




expt" 


1.00 


NA 


-102.7 


Ae+ 


expft 


1.00 


NA 


-126.2 


Ag+/Cl- 


32* 


0.99 


6-pt 


-189.1 


Ag+/Cl- 


expt" 


1.00 


NA 


-191.5 


Ni+ -» Ni2+ 


32 


0.99 


2-pt 


-365.4 


Ni+ ^ Ni^+ 


32 


0.99 


6-pt 


-365.6 


Ni+ ^ Ni^+" 


32 


0.99 


2-pt 


-354.5 


Ni+ ^ Ni^+" 


32 


0.99 


6-pt 


-353.7 



TABLE III: Ag+ hydration free energies, and Ni+ -» Ni^+ 
hydration free energy changes. H2O densities and AGhyd are 
in units of g/cc and kcal/mol, respectively. All simulations 
are based on the PBE functional, except that the DFT-I-U 
formalism with U—i eV is applied for Ni predictions marked 
with an "x." The asterick denotes AGhyd adjusted for packing 
effects. "Ref. [21]. Experimental values adjusted for surface 
potentials are depicted with a dagger; see text for details. 



tial (IP). The VASP PBE PP predicts the Ag IP to be 
178.9 kcal/mol, while the first and second IP for Ni are 
predicted to be 160.6 and 492.9 kcal/mol, respectively. 
Adding the respective AGhyd, Eqs. [5] and [6] yield AG of 
-1-57.5 and -1-76.0 kcal/mol, respectively. These individ- 
ual half-cell reaction AG have not yet been referenced 
to the standard hydrogen potential. The overall Ag+ 
+ Ni+ Ag + Ni^+ reaction, however, does not suf- 
fer from surface potential ambiguities. If we use the IP 
predicted using the PBE functional, the AG of this re- 
action becomes -1-18.5 kcal/mol, or -1-0.80 eV, in water. 
We stress that the pertinent Ag species is the silver atom 
suspended in water, not bulk silver metal. 

PBE predictions for IP are, however, problematic. 
While our pseudopotential PBE method fortuitously pre- 
dicts an Ag IP in reasonable agreement with the experi- 
mental value of 174.6 kcal/mol, the most accurate quan- 
tum chemistry method (CCSD(T)) with relativistic cor- 
rections in fact underestimates this value by ~ 1 eVi^i^ 
While the CCSD(T) method is accurate for the first IP of 
Ni,.^*' our pseudopotential PBE approach severely overes- 
timated the second Ni ionization potential measured at 
418.7 kcal/mol,2i 

A more reasonable approach is to combine experimen- 
tal IP and AIMD AGhyd- This yields AG=+0.01 eV for 
Eq.m The predicted value is significantly more endother- 
mic than the -0.6 eV cited in the experimental radiolysis 
literature i^ii^i^ That -0.6 eV value was derived by esti- 
mating the Ni+ AGhyd using a simple Pauling ionic ra- 
dius and a dielectric continuum approximation;^^ as the 
authors stressed, ligand field effects, which can be a frac- 
tion of an eV for first row transition metal ions in water fS£ 
were neglected. AIMD AGhyd calculations, free from 
these assumptions, should yield more accurate redox po- 



tentials for metal ions in unstable valence states encoun- 
tered as transients in radiolysis experiments i^^'^^i^"^ 

Finally, we note that the Ni^+ AGhyd depends on 
whether the DFT-I-U approach is used in calculating 
{dH{q)/q)q along the AIMD trajectory. Setting [/=4 
(6) cV already decreases the gas phase Ni^+-(Il20)6 clus- 
ter binding energy by ~ 0.5 eV (1.0 eV) without in- 
ducing noticeable changes in the geometry of the com- 
plex. Since the octahedral Ni^+ hydration shell is quite 
stable in liquid water, a similar change in the aqueous 
phase AGhyd is expected if U varies by like amounts. 
We have indeed found that using DFT-f U (C/=4 eV) to 
compute {dH[q)/q)q decreases the solvation by roughly 
12 kcal/mol, yielding a AGhyd of -353.7±1.0 kcal/mol. 
With this DFT-I-U AGhyd, Eq. H becomes endothermic 
by -1-0.51 eV compared with the -1-0.01 eV predicted with 
PBE (i.e., t7=0 eV). Whether PBE or DFT-hU yields 
more accurate AGhyd will be assessed in the future by 
comparison with high level quantum chemistry, new DFT 
functionalSf^ or gas phase experimental values such as 
those reported for monovalent cations and anionsi^ 

The above analysis suggests that predicting re- 
dox potential of half cell electrochemical reactions of 
first row transition metal ions like Ni+ remains a 
challenge r^i^^ and that reported redox values in the radi- 
olysis literature^i^ may need to be extensively revised. 
We stress that our approach, which partitions redox po- 
tentials into hydration free energies and IP, circumvents 
DFT inaccuracies associated with IP predictions. 

IV. CONCLUSIONS 

We have applied ah initio molecular dynamics (AIMD) 
simulations to compute the absolute hydration free ener- 
gies of Li"*" , Cl~ , and Ag+ . While some small contribu- 
tions from packing (entropy) effects and simulation cell 
size dependences for anions still need to be estimated us- 
ing classical force field based simulations, the dominant 
electrostatic contributions come from density functional 
theory (DFT) and rigorous liquid state statistical me- 
chanical methods i^i^iiiiii 

To compare with experimental values, care must 
be taken to account for surface potential con- 
tributions which can be decomposed into water 
dipole and quadrupole ("second spherical moment") 
contributions, ^^'^^ q{4>d + 't'q)- So far, the water- vapor 
interface surface potential has not been computed us- 
ing AIMD/PBE. Nevertheless, the experimental data 
tabulated by Marcus-- and Tissandier et al^ can be 
compared with AIMD values by adding q(j)q and sub- 
tracting q(l)d values estimated using the SPC/E water 
model, respectively. In both cases, we would be com- 
paring with AGhyd values fitted to the SPC/E water 
model; but to the extent that the SPC/E (f)d is an ac- 
curate physical quantity, comparing AIMD AGhyd with 
AGxiss — 0d (SPC/E) (plus a standard state correction 
C*^°^) should be model-independent. With these caveats. 



we find that the AIMD AGhyd for Li+ and C1+ are within 
4.9 (4 %) and 0.4 kcal/mol (0.5 %) of Tissandier et aUs 
values adjusted this way. The deviations from Marcus' 
values, compiled after removing surface potential and 
standard state contributions, are larger, probably due to 
uncertainties in estimates. The sum of AGhyd for the 
Li+/Cl~ ion pair, where surface potential effects cancel, 
agree with the two sets of experimental values to within 
2.3% and 2.6%, respectively;^!^ The Ag+/Cr ion pair 
has a combined AGhyd within 1.2 % of Marcus' data. 

We also compute the change in AGhyd associated with 
Ni"*" being oxidized to Ni^+ . Coupled with the hydration 
free energy of Ag+ and experimental ionization potential 
values, we arrive at a free energy change of 0.01 eV (PBE) 
and 0.51 eV (DFT+U, U=A eV) for the Ag+ + Ni+ ^ 
Ag (atom) + Ni'^^ reaction in water. Whether PBE or 
DFT+U yields more accurate AGhyd will be assessed in 
the future by comparison with high level quantum chem- 
istry, new DFT functionals, or experimental values. This 
calculation is pertinent to predicting the redox poten- 
tial of unstable Ni+ ions. The Ni+ oxidation potential 
often cited in the radiolysis experimental literature actu- 
ally contains a theoretical hydration free energy estimate 
based on the Ni"*" Pauling radius, and it does not account 
for ligand field effects i^i^ Our results suggest that such 
reported values may need to be re-examined with the 
more accurate AIMD approach. 

Even without more accurate determination of surface 
potentials, our formalism can be applied to predict the 
AIMD AGhyd difference between like-charged ions such 
as Na"*" and K"*" , which is relevant to understanding mech- 
anisms of selective ion binding. Our work also paves the 
way for AIMD calculations of the hydration free ener- 
gies of more complex ions and of ions at water-material 
interfaces, inside carbon nanotubes where material po- 
larizability is significant and in inhomogeneous aque- 
ous media in general. Further work on elucidating the 
surface potential entirely with AIMD methods, system- 
atic investigation of the U dependence of hydration free 
energy when DFT-I-U is applied, and comparison with 
other functionals (e.g., BLYP— ) and AIMD packages 
(e.g., CPMD^) will be pursued in the future. 
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FIG. 9: Electrostatic potential as a function of z {(ti{z)) com- 
puted perpendicular to the cubic ice/vacuum interface direc- 
tion and averaged over the lateral {x- and y-) directions. The 
red and black circles depict projections of O and H atoms, 
respectively. 



Appendix: Surface potential from DFT calculations 

[This appendix was not included in the version pub- 
lished in the Journal of Chemical Physics. We thank 
Dr. Shawn Kathmann for useful discussions.] 

A recent DFT calculation2i using the CP2K code^ 
the BLYP exchange-correlation functional, ^'^ and inter- 
facial configurations previously obtained in a 7 ps AIMD 
trajectory^^ has reported a water-vapor surface poten- 
tial of (p ^ —0.018 V. We have not computed (j> directly 
using the VASP— code and the PBE functional. How- 
ever, from the SPC/E and DFT estimates of 4>d and (j)q 
(Ref. lis) respectively, cj) should be -1-4.06 V at 1.0 g/cc 
water density. In the CP2K/BLYP calculation, the wa- 
ter density is ^ 0.9 g/cc, which should yield 90% of (j)q 
at 1.0 g/cc.J^'^^ The functional used in Ref. 74 is also 
different from PBE applied herein. Nevertheless, the 
large difference in (f) values reported is surprising. At 
the same time, (f> depends on the calculation protocol, 
including whether all-electron or pseudopotential meth- 
ods are use d^^i^^ and how long-range electrostatics are 
handled. 

As a preliminary step to a full comparison with Ref. [t^. 
we here use VASP/PBE to directly evaluate <^(z) = 
fx y /LxLy for a cubic ice- vacuum interface. See 

Fig. [9l The ice model contains 64 H2O molecules in a 
12.52x12.52x30.00 simulation ceU. For simplicity, we 
have not relaxed the water geometry and have kept all O- 
H covalent bond at 1.0 A. Since H2O orientational averag- 
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ing is not applicable in crystalline ice, the overall surface 
potential cannot be rigorously decomposed into molec- 
ular contributions . ^^'^^ Also unlike liquid water, trans- 
lational invariance does not apply in a crystal. Thus 
oscillations in (f>{z) emerge as a function of z, with the 
peaks coinciding with the z-planes containing the oxygen 
nuclei. Around the nuclei, V{x, y, z) is at its most pos- 
itive because the nuclear charge in the pseudopotential 
is only weakly screened by valence electrons. (We stress 
that V{x, y, z) is "measured" using a virtual point charge, 
and does not contain Pauli-exclusion repulsion with core 
electrons.) Thus, —\e\V{x,y,z) integrated over x and y 
should be most negative at the nuclear positions. Note 
also that the surface dipole contribution to the overall 
is not negligable. With these caveats, the average dif- 
ference between (j){z) inside the ice crystal and 0(z) in 
the vacuum region is predicted to be -1-2.7 V (Fig. 9), 
which is far closer to the -1-4.06 V we have estimated 
for the PBE liquid/vapor interface at 1.0 g/cc liquid wa- 
ter density than the -0.018 V reported in Ref. UA— We 
have also used VASP/PBE to compute a ~ 3.5 eV for 
a single snapshot of a water-vapor interfacial geometry 
generated using 128 SPC/E H2O molecules (not shown). 
This value is again consistent with the value of (p adopted 
in this work. 

As discussed in Ref. [TqI, 4>q for H2O molecules in liquid 
water is predicted to be 3.85 V using the PBE functional 
and VASP PAW pseudopotentials. It must be on the or- 
der of 4 V, and be positive; if - V, the DFT/PBE 
hydration energy for the Cl~ would be repulsive}^ Look- 
ing at the geometry of the electron cloud of a DFT/PBE 
H2O water molecule versus the point charge distribution 
of the SPC/E model, 0, computed using DFT/PBE must 
indeed carry a sign opposite to that of SPC/E. Given 
the small (j) predicted using the CP2K code^i compared 
to the VASP estimate, it will be interesting to ascer- 
tain whether the intrinsic ion hydration free energies 
AGswaid computed using CP2K are very also different 
from VASP predictions. When added to the CP2K e0, 
CP2K AGEwaid should agree with VASP or experimen- 
tal values. Further effort to reconcile the VASP/PBE 



and CP2K/BLYP surface potentials will be undertaken 
in the future. 

We point out that work functions $ in metals, which 
are the energies required to eject electrons from the ma- 
terials, have been measured and calculated using DFT 
methods, and they typically amount to several electron 
voltsi^ Here $ =e0 — /Xe, where fie is the Fermi level or 
chemical potential of the electron, (j) computed in those 
cases also depend on whether an all-electron or pseudopo- 
tential approach is used, but $ should be independent of 
such details by virtue of cancellation with the /ie term. 

We have so far sidestepped the issue of experimen- 
tal measurement of 0. depends on the salt present 
at the surfacci^ Neither the magnitude nor apparently 
even the sign of this quantity at the water- vapor interface 
has been unambiguously established . ^^1^^ On the other 
hand, to the extent that (j) is mainly of interest to de- 
termining the absolute hydration free energies of ions at 
infinite dilution, it is sufficient to treat as a compu- 
tational method- or force field-dependent entity. How- 
ever, it may be possible to measure 4> directly and un- 
ambiguously in the future. Here we confine ourselves to 
the following comments, (i) Just as the theoretical 
depends on the method used, the measured (j) may be 
sensitive to experimental details, (ii) With optical (e.g., 
sum harmonic generation) measurements where instanta- 
neous electromagnetic fields are involved, the electronic 
structure method used should use as few approximations 
as possible. Thus, all-electron (or frozen core treatment 
of core electrons) calculations and treating the nuclei as 
point charges should yield the most reliable results. Core 
electrons should have minimal impact on most physical 
properties but explicitly add to (jjq-^^ (iii) If an extrapo- 
lation of (j) from the thermodynamic critical point (where 
(j) rigorously vanishes) is applied, it may be useful to note 
that, in electronic structure calculations, the quadrupole 
component 4>q of H2O molecules dominates at moderate 
water density, and that this quantity is proportional to 
Ap, the difference between liquid and vapor densities. 
Thus the slope of d(j)/dT^ where T is the temperature, 
may be most strongly correlated with dAp/dT. 
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Fig. 3 (a): Integrated changes in electron density, A(x) = J dydz[p{x,y, z)n — p{x,y, z)c], as a function of spatial coordinate x 
for the various values of q. p„ and pc axe the densities for the neutral and the charged systems, respectively. All charged 

species, Li''^ have been shifted to a; = 0. Symbols correspond to actual grid-points, the continuous lines are cubic 
interpolations, (b)-(d): Isosurface plots of the electron density difference, p{x,y,z)„ — p{x,y,z)c (iso-value = ± 0.01 a.u., 
white < 0, blue > 0), for q =0.1, 0.6, and 1.0. Periodic boundary conditions apply; the prominent, 8 blue spheres represent 
the (periodically replicated) changes in Li'^ densities, and some changes in water dipole moments are apparent too. See 
Sec. Ill Bi for technical details. (The panels are located at the end of the document.) 



